Final Figures and Stats

PAR

---

Methods

PAR values are from flourecien placed every 5 meters along streams. Three 1mL stint vials were placed every 5 meters in each reach and the average value of the three was used. To account for natural fluctuations in FL flourecence, controls wrapped in tin foil were placed at random 5 meter intervals with the 3 replicate vials. The measured means were then adjusted by the mean natural change in control values.

Figures

FL$Stream <- recode(FL$Stream, Loon = "LOON", MCTE = "MCTE", `W-100` = "W-100", `W-113` = "W-113", Chucksney = "CHUCK")
FL$Stream <- factor(FL$Stream, levels = c("LOON", "MCTE", "W-113", "W-100", "CHUCK"))
FL$Reach<-as.factor(FL$Reach)
FL$Meter<-as.factor(FL$Meter)
FL$PAR<- as.numeric(FL$PAR)

## Warning: NAs introduced by coercion

PAR.mean <- aggregate(PAR ~ Stream + Reach + Meter, mean, data = FL)

PAR.mean %>% filter(Stream != "W-122") %>%
  
  ggplot(aes(x = Meter, y = PAR, color = Reach, group = Reach)) +
    geom_line(size = 3) +
    ylab("PAR (mol/m^2 day^-1)")+ scale_x_discrete(breaks = seq(0, 90, by = 15)) +
    labs(
         y = "PAR (mol m^-2 day^-1)",
         x = "Meter",
         color = "Reach") + 
    ylim(0,18) +
    facet_wrap("Stream") +
    scale_color_viridis_d(option = "E") +
    theme(legend.position = c(.84, .27), legend.key.width = unit(3, "cm"),
          aspect.ratio = 1)

mean(PAR.mean$PAR[PAR.mean$Reach == "Treatment"])

## [1] 3.715202

mean(PAR.mean$PAR[PAR.mean$Reach == "Reference"])

## [1] 1.280928

Chlorophyll \(\alpha\)

---

Methods

Chlorophyll  values were measured using a BenthoTorch. Three tiles every 10 meters

Figures

# Change variables to factors.  Re-order with CHUCK at the end
Both_years$Stream <- factor(Both_years$Stream, levels = c("LOON", "MCTE", "W-113", "W-100", "CHUCK", "W-122"))
Both_years$Year <- factor(Both_years$Year, levels = c("2018", "2017"))
Both_years$Year.Treatment <- as.factor(Both_years$Year.Treatment)
Both_years$Treatment <- as.factor(Both_years$Treatment)
Both_years$Meter <- as.factor(Both_years$Meter)
Both_years$BenthoTotal <- as.numeric(Both_years$BenthoTotal)

# Drop W-122.
Both_years <- Both_years %>% filter(Stream != "W-122")

# Calculate mean of the three tiles for each meter
bentho_both.mean <- aggregate(BenthoTotal ~ Stream + Treatment + Year.Treatment + Meter + Year, mean, data = Both_years)
bentho.reach <- aggregate(BenthoTotal ~ Stream + Treatment + Year.Treatment + Year, mean, data = Both_years)


ggplot(data = bentho_both.mean) +
  geom_line(mapping = aes(x = Meter, y = BenthoTotal, color = Treatment, 
  linetype = Year, group = Year.Treatment), size = 3) +
  labs(
       y = "Chlorophyll a (ug/cm^2)",
       x = "Meter",
       color = "Reach") + 
  scale_x_discrete(breaks = seq(0, 90, by = 15)) +
  scale_color_viridis_d(option = "E") +
  facet_wrap("Stream") +
  theme(legend.position = c(.86, .2), legend.key.width = unit(3, "cm"),
        aspect.ratio = 1)

Benthic Invertebrates

---

Benthic Setup

---

In order to calculate density we must take into account that our count values are from a subsample of three pooled samples.

#Calculations...
bugs$Density <- bugs$Count / bugs$PercentSub / .09  #Calculate total in pooled sample

bugs$CollDate <- as.Date(bugs$CollDate, format = "%m/%d/%y") %>% year() %>% as.factor()

bugs.agg <- bugs %>%
  dplyr::group_by(CollDate, Stream, Treatment, Taxon) %>%
  summarise_at(vars(Density), sum) %>% ungroup

bugs.agg$Density <- bugs.agg$Density / 3    #Divide by number of samples pooled

---

We divide Count by percent subsampled (actually a fraction) to get a count for the total sample taken. We then divide by .09 which is the area of the surber sampler (in m\(^2\)). We then aggregate and divide by three (There were three samples taken per reach).

# Spread and then gather to fill missing Taxons with 0's
bugs.agg <- spread(bugs.agg, key = "Taxon", value = "Density") %>% #bugs.agg is benthic samples
  mutate_if(is.numeric , replace_na, replace = 0) %>% 
  gather(key = "Taxon", value = "Density", 4:ncol(.))

# Add FFG and size classes
bugs.agg <- merge(bugs.agg, sizes) %>% merge(., ffg)

bugs.agg <- bugs.agg %>% mutate(StreamTreat = interaction(Stream, Treatment, sep = " "))
bugs.agg <- bugs.agg %>% mutate(YearTreat = interaction(CollDate, Treatment, sep = " "))

# Create list of CollDate, Stream and Taxon
bugs.reach.diff <- bugs.agg %>% select(-c("Density", "CollDate")) %>% unique()

# Calculate Differences
bugs.reach.diff$Diff <- bugs.agg$Density[bugs.agg$CollDate == "2018"] -     
  bugs.agg$Density[bugs.agg$CollDate == "2017"]

# Filter for only 2018 bugs
bugs18 <- filter(bugs.agg, CollDate == "2018")

bugs.treat <- bugs.agg %>% 
  group_by(Stream, Treatment, StreamTreat, YearTreat, FFG) %>%
  summarise_at(vars(Density), funs(sum)) %>% ungroup

## Warning: funs() is soft deprecated as of dplyr 0.8.0
## please use list() instead
## 
##   # Before:
##   funs(name = f(.))
## 
##   # After: 
##   list(name = ~ f(.))
## This warning is displayed once per session.

bugs.treat <- bugs.agg %>% select(-c("Density", "CollDate")) %>% unique()

Statistics

# Find total density for all bugs per stream-reach and total density by FFG
bug_totals <- bugs.agg %>% dplyr::group_by(Stream, Treatment, CollDate, YearTreat, StreamTreat) %>% 
  summarise_at(vars(Density), sum)
ffg_totals <- bugs.agg %>% dplyr::group_by(Stream, Treatment, CollDate, YearTreat, StreamTreat, FFG) %>% 
  summarise_at(vars(Density), sum)

# Log transform density
bug_totals <- bug_totals %>% mutate(logDensity = log(Density))
ffg_totals <- ffg_totals %>% mutate(logDensity = log(Density))

We will start with plotting total invertebrate density for a single stream with both reaches for both years in order to get a feel for the raw data. Ultimately we will be taking the ratio between years for each reach and then comparing the change between years in the treatment reach to the change between years of the control reach to assess the effects of the canopy-opening treatment.

bug_totals %>% filter(Stream == "W-113") %>%
  ggplot(aes(x = YearTreat, y = Density, group = CollDate, fill = Treatment)) +
  geom_col() +
  labs(title = "W-113 Invertebrate Density by Year and Reach",
       x = "Year and Treatment \n (Y = Treatment reach, N = Control reach)",
       y = "Invertebrate Density") +
    scale_fill_viridis_d(option = "E", end = .9)

Next we plot log ratios between years (2018Y / 2017Y and 2018N / 2017N) for a single stream

A log ratio of 0 indicates no change in abundance between 2017 and 2018. Positive values mean an increase in abundance in the post-treatment year, and a negative value would suggest a decline in abundance in the post-treatment year. We then compare the yearly ratio’s of each reach (control or treatment) to assess the effect of the treatment.

ffg.ratios$Year <- as.factor(ffg.ratios$Year)
ffg.ratios %>%  group_by(Stream, Year) %>% 
  summarise_at(vars(Difference, Ratio, `Log(Ratio) Year`), mean) %>% 
  
  filter(Stream == "W-113") %>% 
  ggplot(aes(x = Year, y = Difference, fill = Year)) +
  geom_col() +
  labs(title = "W-113 Invertebrate Density Differences Between Reaches",
       x = "Year",
       y = "Difference of Invert Density (Treatment - Control)") +
  scale_fill_viridis_d(option = "E", end = .9) +
  geom_hline(aes(yintercept = 0))

Now we plot all the log ratios

ffg.ratios %>%  group_by(Stream, Year) %>% 
  summarise_at(vars(Difference), mean) %>% 
  spread(key = Year, value = Difference) %>% mutate(ratio = `2018` / `2017`) %>%
  
  ggplot(aes(y = ratio, x = Stream)) +
  geom_col() +
  labs(title = "W-113 Log of Invertebrate Density Ratio's \n (2018 / 2017)",
       x = "Year",
       y = "Reach Difference of Invert Density") +
  scale_fill_grey() +
    scale_color_grey()

And now the average log ratio across streams for each reach with 95% confidence interval calculated using the t score from Welch’s t-test.

bug.ratios <- ffg.ratios %>% group_by(Stream, Year) %>% 
  summarise_at(vars(Difference, Ratio, `Log(Ratio) Year`), mean)

bug.ratios$Year <- as.factor(bug.ratios$Year)

#t.test(Difference ~ Year, data = bug.ratios)

avg.ratio <- ddply(bug.ratios, .(Year), summarize,
                   Mean = mean(`Log(Ratio) Year`),
                   SD = sd(`Log(Ratio) Year`))

avg.difference <- ddply(bug.ratios, .(Year), summarize,
                   Mean = mean(`Difference`),
                   SD = sd(`Difference`))


avg.difference %>% 
  ggplot(aes(x = Year, y = as.numeric(Mean), color = Year)) +
  geom_point(position = "dodge") +
  geom_errorbar(aes(x = Year, ymin = Mean - ((SD*2.345)/sqrt(5)), ymax = Mean + ((SD*2.345)/sqrt(5)), #2.345 is the critical value 7.2983 degrees of freedom 
                    color = Year), position = "dodge") +
  labs(title = "Average of the Reach Difference",
       x = "Year",
       y = "Difference") +
  theme_bw(base_family = "LM Roman 10") +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold")) 

## Warning: Width not defined. Set with `position_dodge(width = ?)`

And now do it all again for FFG’s. We start with looking at one stream, then the log ratio for that stream, then the log ratio for all streams and finally the average log ratio for each FFG with error bars of one SD.

ffg.ratios %>%  group_by(Stream, Year) %>% 
  filter(Stream == "W-113") %>% 
  ggplot(aes(x = Year, y = Difference, fill = Year)) +
  geom_col() +
  labs(title = "W-113 Log of Invertebrate FFG Density Ratio's \n (2018 / 2017)",
       x = "Year",
       y = "Difference") +
  facet_wrap("FFG") +
  scale_fill_grey() +
    scale_color_grey()

ffg.ratios %>%  
  ggplot(aes(x = Stream, y = Difference, fill = Year)) +
  geom_col(position = "dodge") +
  labs(title = "Invert Density Difference",
       x = "Year",
       y = "Difference") +
  theme_bw(base_family = "LM Roman 10") +
  theme(legend.position = c(.1, .4), 
          legend.key.width = unit(5, "cm"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold")) +
  facet_wrap("FFG") +
  scale_fill_grey()+
    scale_color_grey()

ffg.avg.ratio <- ddply(ffg.ratios, .(FFG, Year), summarize,
      Mean = mean(`Log(Ratio) Year`),
      SD = sd(`Log(Ratio) Year`))

ffg.avg.diff <- ddply(ffg.ratios, .(FFG, Year), summarize,
      Mean = mean(`Difference`),
      SD = sd(`Difference`))

ffg.avg.diff$Year <- as.factor(ffg.avg.diff$Year)

ffg_names <- c("SH", "P", "SCe", "CG", "SCi", "CF")

ffg_results <- data.frame("FFG" = numeric(6), "tval" = numeric(6), "pval" = numeric(6))  #initialize data frame
for (i in ffg_names){
  k = match(i, ffg_names)  #find the index of i in ffg_names
  x <- filter(ffg.ratios, FFG == i) %>%
  t.test(`Difference` ~ Year, data = .)
  ffg_results[k, "FFG"] <- i
  ffg_results[k, "tval"] <- x$statistic
  ffg_results[k, "pval"] <- x$p.value
}


ffg.avg.diff %>% 
  ggplot(aes(x = FFG, y = Mean, fill = Year)) +
  geom_col(position = "dodge") +
  geom_errorbar(aes(x = FFG, ymin = Mean - ((SD*2.345)/sqrt(5)), ymax = Mean + ((SD*2.345)/sqrt(5)), 
                    color = Year), width = .7, 
                    position = position_dodge(.9)) +
  labs(
       x = "Functional Feeding Group",
       y = bquote('Density Difference'~("Treatment - Control; "~m^-2))) +
 theme_bw(base_family = "LM Roman 10") +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold")) +
        geom_hline(yintercept = 0, color = "black", size = .4) +
  scale_fill_viridis_d(option = "E", end = .9) +
    scale_color_viridis_d(option = "E", end = .9) 

NMS

---

Who knows what I’m actually doing, all of this ordination stuff is currently pre BOT 570. I will definitely update all this once I know how to handle singleton taxa, loads of zero’s, etc.

# Spread data to get a taxon matrix, create factor of CollDate and Treatment
bugs.mds <- bugs.agg %>% 
  select(-c("FFG", "Size", "EPT")) %>%
  group_by(Stream, Treatment, CollDate, Taxon) %>%  
  spread(key = "Taxon", value = "Density") %>%
  mutate_if(is.numeric , replace_na, replace = 0) %>% ungroup() %>%
  mutate(YearTreat = interaction(CollDate, Treatment) %>% factor(levels = c("2017.N", "2017.Y", "2018.N", "2018.Y"))) 

## `mutate_if()` ignored the following grouping variables:
## Columns `Stream`, `Treatment`, `CollDate`

# Get grouping variables and taxon matrix
Density <- bugs.mds %>% 
  select("Ameletus":"Zapada")

env_bugs <- Density %>% select_if(funs(sum(.) > 200))

grouping <- select(Enviro, Stream, Treatment, CollDate, YearTreat, YearTreatQ, Chla, PAR) %>% cbind(env_bugs)

## Run 0 stress 0.09528951 
## Run 1 stress 0.1049411 
## Run 2 stress 0.09528901 
## ... New best solution
## ... Procrustes: rmse 0.0002408672  max resid 0.0005452333 
## ... Similar to previous best
## Run 3 stress 0.09528869 
## ... New best solution
## ... Procrustes: rmse 0.0001598907  max resid 0.0003153504 
## ... Similar to previous best
## Run 4 stress 0.1074884 
## Run 5 stress 0.0980021 
## Run 6 stress 0.0952887 
## ... Procrustes: rmse 3.875998e-05  max resid 8.33822e-05 
## ... Similar to previous best
## Run 7 stress 0.1006332 
## Run 8 stress 0.09448684 
## ... New best solution
## ... Procrustes: rmse 0.03088568  max resid 0.1057163 
## Run 9 stress 0.09746273 
## Run 10 stress 0.09894043 
## Run 11 stress 0.09746193 
## Run 12 stress 0.09448911 
## ... Procrustes: rmse 0.0009165532  max resid 0.003273089 
## ... Similar to previous best
## Run 13 stress 0.104941 
## Run 14 stress 0.09895678 
## Run 15 stress 0.09528933 
## Run 16 stress 0.09746181 
## Run 17 stress 0.0944871 
## ... Procrustes: rmse 0.0004116904  max resid 0.001470199 
## ... Similar to previous best
## Run 18 stress 0.09895678 
## Run 19 stress 0.09895399 
## Run 20 stress 0.0989412 
## *** Solution reached

#set up NMDS with dimensions of sol and env factors from "Enviro". 
NMDS <- data.frame(x = sol$points[, 1], y = sol$points[ ,2], 
                   Stream = select(grouping, Stream), 
                   Treatment = select(grouping, Treatment), 
                   CollDate = select(grouping, CollDate),
                   YearTreat = select(grouping, YearTreat),
                   Chla = select(grouping, Chla),
                   YearTreatQ = select(grouping, YearTreatQ),
                   PAR = select(grouping, PAR))

Reach17N <- NMDS[NMDS$YearTreat == "2017.N", ][chull(NMDS[NMDS$YearTreat == 
    "2017.N", c("x", "y")]), ] 

Reach17Y <- NMDS[NMDS$YearTreat == "2017.Y", ][chull(NMDS[NMDS$YearTreat == 
    "2017.Y", c("x", "y")]), ] 

Reach18N <- NMDS[NMDS$YearTreat == "2018.N", ][chull(NMDS[NMDS$YearTreat == 
    "2018.N", c("x", "y")]), ] 

Reach18Y <- NMDS[NMDS$YearTreat == "2018.Y", ][chull(NMDS[NMDS$YearTreat == 
    "2018.Y", c("x", "y")]), ] 

hull.data <- rbind(Reach17N, Reach17Y, Reach18Y, Reach18N) 

vec.env <- envfit(sol, grouping, perm = 1000, na.rm = TRUE)
vec.env.df <- scores(vec.env, display = "vectors") %>% as_tibble()
vec.env.df$r <- vec.env$vectors$r
vec.env.df$pval <- vec.env$vectors$pvals
vec.env.df$names <- colnames(grouping[-c(1,2,4)])
vec.env.df <- vec.env.df %>% filter(r > .36) %>% slice(2:3)

ggplot(data = NMDS, aes(x, y)) +
  geom_polygon(hull.data, mapping = aes(x = x, y = y, color = YearTreat, group = YearTreat), size = 2, alpha = 0) +
  geom_point(aes(x = x, y = y, shape = Stream), size = 5) +
  geom_segment(data = vec.env.df, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
      arrow = arrow(length = unit(0, "cm")), colour = "red", inherit_aes = FALSE, size = 1) +
  geom_text(data = vec.env.df, aes(x = NMDS1, y = NMDS2, label = names), size = 8, color = "red") +
  geom_text_repel(data = NMDS, aes(x = x, y = y, label = Stream), size = 6) +
  labs(
       x = "Axis 1",
       y = "Axis 2") +
  scale_color_viridis_d() +
  theme_bw(base_family = "LM Roman 10") +
    theme(axis.text.x = element_blank(),  # remove x-axis text
        axis.text.y = element_blank(), # remove y-axis text
        axis.ticks = element_blank(),  # remove axis ticks
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank(),
        legend.position = c(.08, .76)) 

## Warning: Ignoring unknown parameters: inherit_aes

ed_scraper_totals <- ffg_totals %>% filter(FFG == "SCe")
in_scraper_totals <- ffg_totals %>% filter(FFG == "SCi")
scraper_totals <- rbind(ed_scraper_totals, in_scraper_totals)
scraper_agg <- scraper_totals %>% group_by(Stream, Treatment, YearTreat, CollDate) %>% summarise_at(vars(Density, logDensity), .funs = sum)

Actually maybe cool result, all of the treatment reaches now plot closer post-gap year. The controls show a similar pattern though maybe a little less pronounced. Everything also shifted left between the two years, although this isn’t because of the treatment, just annual changes. Maybe the gap amplified whatever effect the year had, extra, extra light or temperature maybe.

Summary Info

---

bugs.reach.diff %>%
group_by(Stream, Treatment, FFG) %>%
  summarise_at(vars(Diff), funs(sum)) %>%
  
  ggplot(aes(x = FFG, y = Diff, fill = Treatment)) + 
  geom_col(position = "dodge") +
  ggtitle("Total Inverebrate Density (m^2)") +
  scale_fill_viridis_d(option = "E") +
  facet_wrap("Stream")

The gap has no real effect, except maybe MCTE. We saw this in the t tests and bar plots, and in the ordination.

Diets

---

Diet Setup

---

# Fill missing family's with order
Diets$Family[is.na(Diets$Family)] <- Diets$Order[is.na(Diets$Family)]

# Merge diets, sizes and FFG's
Diets <- merge(Diets, sizes, by.x = "Family", by.y = "Taxon") %>% 
  merge(., ffg, by.x = "Family", by.y = "Taxon")

# Aggregate by rep
Diets.by.fish <- Diets %>% 
  group_by(Family, Stream, Treatment, Origin, FFG, Size, Rep) %>%
  summarise_at("Count", funs(sum)) %>% ungroup

# Aggregate by reach
Diets.reach <- Diets %>% 
  group_by(Family, Stream, Treatment, Origin, FFG, Size) %>%
  summarise_at("Count", funs(sum)) %>% ungroup

# I, Ivlev's Selectivity Index
bugndiet.ffg <- bugndiet.ffg %>% mutate(I = (frac.diet - frac.benthic) / (frac.diet + frac.benthic))

# L, Linear selectivity index
bugndiet.ffg <- bugndiet.ffg %>% mutate(L = (frac.diet - frac.benthic))

# D, Difference selectivity index
bugndiet.ffg <- bugndiet.ffg %>% mutate(D = 
                          (frac.diet - frac.benthic) / (frac.diet + frac.benthic - 2*frac.diet*frac.benthic))

# Q, Quotient selectivity index
bugndiet.ffg <- bugndiet.ffg %>% mutate(Q = (frac.diet*(1-frac.benthic)) / (frac.benthic*(1-frac.diet)))
      

# Calculate average proportion of FFG's in benthos and diet of all streams and avg index values
avg.ffg.diet <- ddply(bugndiet.ffg, .(Treatment, FFG), summarize,
                   Mean.benthic = mean(frac.benthic),
                   SD.benthic = sd(frac.benthic),
                   Mean.diet = mean(frac.diet),
                   SD.diet = sd(frac.diet),
                   Mean.I = mean(I),
                   SD.I = sd(I),
                   Mean.L = mean(L),
                   SD.L = sd(L),
                    Mean.D = mean(D),
                   SD.D = sd(D),
                   Mean.Q = mean(log(Q)),
                   SD.Q = sd(log(Q)))


# Select only "important" families
bugndiet.select.fam <- bugndiet.fam %>% mutate(Family = recode(Family, "Brachycentridae" = "Brachy", "Chironomidae" = "Chiro", "Juga" = "Juga", "Baetidae" = "Baetid", "Perlidae" = "Perlid", "Ostrocod" = "Ostro", "Elmidae" = "Elmid", "Rhyacophilidae" = "Rhyaco", "Dixidae" = "Dixid", .default = "")) %>% filter(Family != "")

# I, Ivlev's Selectivity index
bugndiet.select.fam <- bugndiet.select.fam %>% mutate(I = (frac.diet - frac.benthic) / (frac.diet + frac.benthic))

# L, Linear selection Index
bugndiet.select.fam <- bugndiet.select.fam %>% mutate(L = (frac.diet - frac.benthic))

# D, Difference selectivity index
bugndiet.select.fam <- bugndiet.select.fam %>% mutate(D = 
                          (frac.diet - frac.benthic) / (frac.diet + frac.benthic - 2*frac.diet*frac.benthic))

# Q, Quotient selectivity index
bugndiet.select.fam <- bugndiet.select.fam %>% mutate(Q = (frac.diet*(1-frac.benthic)) / (frac.benthic*(1-frac.diet)))
                      

# Calculate average proportion of Families in benthos and diet of all streams
avg.fam.diet <- ddply(bugndiet.select.fam, .(Treatment, Family), summarize,
                   Mean.benthic = mean(frac.benthic),
                   SD.benthic = sd(frac.benthic),
                   Mean.diet = mean(frac.diet),
                   SD.diet = sd(frac.diet),
                   Mean.I = mean(I),
                   SD.I = sd(I),
                   Mean.L = mean(L),
                   SD.L = sd(L),
                   Mean.D = mean(D),
                   SD.D = sd(D),
                   Mean.Q = mean(log(Q)),
                   SD.Q = sd(log(Q)))

Diets.reach$Stream.Treat <- interaction(Diets.reach$Treatment, Diets.reach$Stream)

Diets.reach %>%
  group_by(Stream.Treat, Origin) %>%
  summarise_at(vars(Count), funs(sum)) %>% ungroup() %>%

ggplot(aes(x = Stream.Treat, y = Count, fill = Origin)) +
  geom_col(position = "stack") +
  labs(title = "Mean Aquatics and Terrestrials \n in Diets by Reach",
       y = "Mean # in Diet",
       x = "Stream and Reach") +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    scale_fill_viridis_d(option = "E", end = .8)

By Family

---

ggplot(aes(x = Family, y = Count, fill = Treatment), data = Diets.reach) +
  geom_col() +
  coord_flip() 

fam.agg <- Diets.reach %>% group_by(Family) %>% summarise_at(vars(Count), sum)
fam.agg$Count[fam.agg$Family == "Chironomidae"] / sum(fam.agg$Count)

## [1] 0.1837161

# Plot proportion of a taxon in the benthic community versus proporion of a taxon in the diet community
ggplot(aes(x = frac.benthic, y = frac.diet, color = Treatment), data = bugndiet.select.fam) +
    geom_point(size = 8) +
    geom_abline(slope = 1, intercept = 0) +
    coord_cartesian(ylim = c(0, .8), xlim = c(0, .8)) +
    geom_text_repel(aes(x = frac.benthic, y = frac.diet, label = Family), size = 8, point.padding = 1, box.padding = .5) +
    labs(title = "Prop. of Taxa in Benthic Comm. \n Vs. Prop. in Diets",
         x = "Proportion of Benthic",
         y = "Proportion of Diet") +
    facet_wrap("Stream", ncol = 3) +
    scale_color_viridis_d(option = "E") +
    theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_rect(fill = "white", color = "grey"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold"),
          legend.position = c(.84, .27), legend.key.width = unit(3, "cm"),
          aspect.ratio = 1)

ggplot(aes(x = Mean.benthic, y = Mean.diet), data = avg.fam.diet) +
    geom_jitter(aes(color = Treatment)) +
    geom_abline(slope = 1, intercept = 0, color = "grey") +
    coord_cartesian(ylim = c(0, .47), xlim = c(0, .47)) +
    geom_text_repel(aes(x = Mean.benthic, y = Mean.diet, label = Family), 
                    force = 1.5, segment.alpha = .5, segment.color = "grey") +
    labs(title = "Prop. of Families in Benthic Comm. Vs. Prop. in Diets", 
         x = "Proportion of Benthic",
         y = "Proportion of Diet") +
    scale_color_viridis_d(option = "E", end = .8) +
    theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_rect(fill = "white", color = "grey"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold"),
          aspect.ratio = 1,
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

ggplot(aes(x = Family, y = Mean.D, color = Treatment, group = Treatment), data = avg.fam.diet) +
  geom_point(position = position_dodge(width = .6)) +
  geom_errorbar(aes(ymin = Mean.D - SD.D, ymax = Mean.D + SD.D, color = Treatment), position = position_dodge(width = .6)) +
    labs( 
         x = "Family",
         y = "Difference Selection Index") +
    scale_color_viridis_d(option = "E", end = .9) +
  theme_bw(base_family = "LM Roman 10") +
    theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_rect(fill = "white", color = "grey"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold"),
          aspect.ratio = 1,
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

Plot of proportional abundance of an individual taxon in the benthic community versus the proportional abundance of the aggregate diets for that reach.

Here we see that Chironomidae typically make up the greatest abundance of both the fish diets and the benthic community, with LOON being the exception (Brachycentridae are equally represented in the diets of control reach and over represented in the treatment). In MCTE, Juga compose a large part of the diets relative to their abundance due to one outlier fish.

There are no overarching differences in how a taxa maps out given the treatment, but the difference in Brachycentridae in LOON between treatment and control is interesting.

By FFG

---

# Plot proportion of a FFG in the benthic community versus proporion of a FFG in the diet community
bugndiet.ffg %>% 

  ggplot(aes(x = frac.benthic, y = frac.diet, color = Treatment)) +
    geom_point() +
    geom_abline(slope = 1, intercept = 0) +
    coord_cartesian(ylim = c(0, 1), xlim = c(0, 1)) +
    geom_text_repel(aes(x = frac.benthic, y = frac.diet, label = FFG), point.padding = .5, box.padding = .5) +
    facet_wrap("Stream", ncol = 4) +
    labs(title = "Prop. of FFG's in Benthic Comm. Vs. Prop. in Diets", 
         x = "Proportion of Benthic",
         y = "Proportion of Diet") +
    scale_color_viridis_d(option = "E") +
    theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_rect(fill = "white", color = "grey"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold"),
        legend.position = c(.8, .15),
          axis.title.x = element_text(hjust = -.1),
          aspect.ratio = 1,
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

Plot of proportional abundance of an individual FFG in the benthic community versus the proportional abundance of the aggregate diets for that reach.

We see that Collector Gatherers and Shredders are the most abundant both in diets and in the benthic community. There aren’t any over arching trends in how FFG’s fall out by treatment.

ggplot(aes(x = Mean.benthic, y = Mean.diet, color = Treatment), data = avg.ffg.diet) +
    geom_point() +
    geom_abline(slope = 1, intercept = 0, color = "grey") +
    coord_cartesian(ylim = c(0, .6), xlim = c(0, .6)) +
    geom_text_repel(aes(x = Mean.benthic, y = Mean.diet, label = FFG), 
                    force = 1.5, segment.alpha = .5, segment.color = "grey") +
    labs(title = "Prop. of FFG's in Benthic Comm. Vs. Prop. in Diets", 
         x = "Proportion of Benthic",
         y = "Proportion of Diet") +
    scale_color_viridis_d(option = "E", end = .8) +
    theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_rect(fill = "white", color = "grey"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold"),
          aspect.ratio = 1,
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

ggplot(aes(x = FFG, y = Mean.L, color = Treatment), data = avg.ffg.diet) +
  geom_point(position = position_dodge(width = .6)) +
  geom_errorbar(aes(ymin = Mean.L - SD.L, ymax = Mean.L + SD.L, color = Treatment), position = position_dodge(width = .6)) +
    labs( 
         x = "Functional Feeding Group",
         y = "Linear Selection Index") +
    scale_color_viridis_d(option = "E", end = .9) +
  theme_bw(base_family = "LM Roman 10") +
    theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_rect(fill = "white", color = "grey"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold"),
          aspect.ratio = 1,
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

ggplot(aes(x = FFG, y = Mean.D, color = Treatment), data = avg.ffg.diet) +
  geom_point(position = position_dodge(width = .6)) +
  geom_errorbar(aes(ymin = Mean.D - SD.D, ymax = Mean.D + SD.D, color = Treatment), position = position_dodge(width = .6)) +
    labs(
         x = "Functional Feeding Group",
         y = "Difference Selection Index") +
    scale_color_viridis_d(option = "E", end = .9) +
  theme_bw(base_family = "LM Roman 10") +
    theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_rect(fill = "white", color = "grey"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
        legend.title = element_text(face = "bold"),
          aspect.ratio = 1,
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

Costello method plot

---

# Find fraction of each family in each individual stomach 
#diets.fish <- ddply(Diets.by.fish, .(Treatment, Stream, Rep), summarize,
                  # frac.diet = Count / sum(Count),
                 # Family = Family)

#family_names <- diets.fish$Family

# Still working this out

#diets.fish <- ddply(Diets.by.fish, .(Treatment, Stream, Rep), summarize,
                  
                  #Family = Family)

#frac.diet[k] <- for(i in family_names){
                    #length(diets.fish$Rep[diets.fish$Family == i]) / max(diets.fish$Rep) 
                  # }

# Calculate the number of reps divided by the total number of reps a family appears in for each reach
#Diets.fish <- Diets.fish %>% 
  #group_by(Stream, Treatment, Family) %>%
  #mutate(frac.fish = length(Family[Count > 0])) %>%
  #group_by(Stream, Treatment) %>%
 # mutate(frac.fish = frac.fish / max(Rep)) %>%
 # group_by(Stream, Treatment, Family)

#Diets.fish <- Diets.fish %>%
  #summarise_at(vars(frac.fish, frac.diet), funs(mean))

#ggplot(data = Diets.fish) +
 # geom_point(mapping = aes(x = frac.fish, y = frac.diet, color = Treatment)) +
 # geom_text(aes(x = frac.fish, y = frac.diet, label = Family), size = 3, nudge_y = .02) +
 # geom_abline(slope = 1, intercept = 0) +
 # coord_cartesian(ylim = c(0, 1), xlim = c(0, 1)) +
 # facet_wrap("Stream") +
  #ggtitle("Costello Plot of Prop. Fish Occurance Vs. Prop. Diet Community") +
 # scale_color_viridis_d(option = "E") +
 # theme(legend.position = c(.84, .27), legend.key.width = unit(5, "cm"))

The Costello method plots % occurance of a taxon in fish versus % of aggregate diet. With the adjusted method the aggregate diet is only of fish that had the taxon present (ignoring zero values). We see that in MCTE and in W-100 there seem to be just a couple fish that specifically target Juga and ELmidae respectively. In MCTE, the taxa from treatment diets seem to consistently constitute a smaller portion of the diet than would be expected given how many fish they occur in.

Diet MDS

---

Who knows what I’m actually doing, all of this ordination stuff is currently pre BOT 570. I will definitely update all this once I know how to handle singleton taxa, loads of zero’s, etc.

# Spread data to get a taxon matrix, create factor of CollDate and Treatment
bugs.mds <- Diets.reach %>% 
  filter(Origin %in% c("A", "T")) %>%
  select(-c("FFG", "Size", "Origin")) %>%
  slice(-65) %>%
  spread(key = "Family", value = "Count") %>%
  mutate_if(is.numeric , replace_na, replace = 0) 

# Get grouping variables and taxon matrix
Enviro <- bugs.mds %>%
  select(Stream, Treatment) 

Density <- bugs.mds %>% 
  select("Amphizoidae":"Uenoidae")

grouping <- select(Enviro, Stream, Treatment)

## Wisconsin double standardization
## Run 0 stress 0.128467 
## Run 1 stress 0.1284675 
## ... Procrustes: rmse 0.0002980788  max resid 0.0007118749 
## ... Similar to previous best
## Run 2 stress 0.1284669 
## ... New best solution
## ... Procrustes: rmse 0.0001467207  max resid 0.000308508 
## ... Similar to previous best
## Run 3 stress 0.1204532 
## ... New best solution
## ... Procrustes: rmse 0.2198124  max resid 0.4822997 
## Run 4 stress 0.1204536 
## ... Procrustes: rmse 7.893107e-05  max resid 0.0001206914 
## ... Similar to previous best
## Run 5 stress 0.1284668 
## Run 6 stress 0.1204532 
## ... New best solution
## ... Procrustes: rmse 0.0001527311  max resid 0.0003160475 
## ... Similar to previous best
## Run 7 stress 0.1204539 
## ... Procrustes: rmse 0.0005409296  max resid 0.001147682 
## ... Similar to previous best
## Run 8 stress 0.3314815 
## Run 9 stress 0.1204531 
## ... New best solution
## ... Procrustes: rmse 6.453872e-05  max resid 0.0001474859 
## ... Similar to previous best
## Run 10 stress 0.1204536 
## ... Procrustes: rmse 0.0002944357  max resid 0.0006053784 
## ... Similar to previous best
## Run 11 stress 0.1204532 
## ... Procrustes: rmse 0.0001322216  max resid 0.0002723085 
## ... Similar to previous best
## Run 12 stress 0.1204532 
## ... Procrustes: rmse 5.347416e-05  max resid 0.0001046169 
## ... Similar to previous best
## Run 13 stress 0.1284666 
## Run 14 stress 0.3195667 
## Run 15 stress 0.1204532 
## ... Procrustes: rmse 9.38903e-05  max resid 0.0001466763 
## ... Similar to previous best
## Run 16 stress 0.1284667 
## Run 17 stress 0.1284667 
## Run 18 stress 0.2194268 
## Run 19 stress 0.1284669 
## Run 20 stress 0.1204534 
## ... Procrustes: rmse 0.0002883572  max resid 0.0005966145 
## ... Similar to previous best
## *** Solution reached

#set up NMDS with dimensions of sol and env factors from "Enviro". 
NMDS <- data.frame(x = sol$points[, 1], y = sol$points[ ,2], 
                   Stream = select(grouping, Stream), 
                   Treatment = select(grouping, Treatment))

Reach18N <- NMDS[NMDS$Treatment == "2018.N", ][chull(NMDS[NMDS$Treatment == 
    "2018.N", c("x", "y")]), ] 

Reach18Y <- NMDS[NMDS$Treatment == "2018.Y", ][chull(NMDS[NMDS$Treatment == 
    "2018.Y", c("x", "y")]), ] 

hull.data <- rbind(Reach18Y, Reach18N) 

ggplot(data = NMDS, aes(x, y, color = Treatment)) +
  geom_polygon(hull.data, mapping = aes(x = x, y = y, fill = Treatment, group = Treatment), alpha = 0) +
  annotate("text", x = (NMDS$x), y = (NMDS$y) + .04, label = NMDS$Stream, size = 3) +
  geom_point(aes(shape = Stream)) +
  ggtitle("MDS of Community in Fish Diets") +
  scale_color_viridis_d() +
  theme(axis.text.x = element_blank(),  # remove x-axis text
        axis.text.y = element_blank(), # remove y-axis text
        axis.ticks = element_blank(),  # remove axis ticks
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank())

  #theme(legend.position = c(.84, .27), legend.key.width = unit(5, "cm"))

There seems to be total overlap of the diet communities in the control versus treatment reaches. This seems to be consistent with what we saw in other plots.

diet.table <- merge(Diets, ffg, by.x = "Family", by.y = "Taxon")

ept <- bugs.agg
ept <- ept %>% filter(Density > 0 ) %>% dplyr::group_by(Stream, Treatment, CollDate) %>% summarise_at(vars(EPT), funs(sum), na.rm = T) %>% ungroup()
ept <- ept %>% spread(Treatment, EPT)
ept$rat <- ept$`Y` / ept$`N`
ept <- ept %>% select(-c(`Y`, `N`))
ept <- ept %>% spread(CollDate, rat)
ept$ddiff <- ept$`2018` / ept$`2017`

Snail <- merge(snails, bentho_both.mean, by.x = c("Stream", "Year", "Treatment", "Meter"), by.y = c("Stream", "Year", "Treatment", "Meter"))

ggplot(aes(BenthoTotal, Snail), data = Snail) +
  geom_point() +
  geom_smooth(method = "lm")

lm.snail <- lm(BenthoTotal ~ Snail, data = Snail)
summary(lm.snail)

## 
## Call:
## lm(formula = BenthoTotal ~ Snail, data = Snail)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.24829 -0.13942 -0.04114  0.10734  0.53289 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.236141   0.023389  10.096 8.24e-16 ***
## Snail       0.032742   0.005538   5.912 8.43e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1782 on 78 degrees of freedom
## Multiple R-squared:  0.3095, Adjusted R-squared:  0.3006 
## F-statistic: 34.96 on 1 and 78 DF,  p-value: 8.432e-08

stream_snail <- snails %>% group_by(Stream, Treatment, Year) %>%
            summarise_at(vars(Snail), sum)

snail_baci <- stream_snail %>% spread(key = Treatment, value = Snail) %>% mutate(diff = Y - N) %>% select(-c("N", "Y"))

Citations

---

Session Info:

sessionInfo()

## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.14.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] plyr_1.8.4        here_0.1          ggrepel_0.8.0    
##  [4] DT_0.5            viridis_0.5.1     viridisLite_0.3.0
##  [7] ggthemes_4.0.1    readxl_1.1.0      lubridate_1.7.4  
## [10] vegan_2.5-2       lattice_0.20-35   permute_0.9-4    
## [13] forcats_0.3.0     stringr_1.4.0     dplyr_0.8.1      
## [16] purrr_0.3.2       readr_1.1.1       tidyr_0.8.1      
## [19] tibble_2.1.1      ggplot2_3.1.1     tidyverse_1.2.1  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.1       rprojroot_1.3-2  assertthat_0.2.1 digest_0.6.19   
##  [5] R6_2.4.0         cellranger_1.1.0 backports_1.1.2  evaluate_0.13   
##  [9] httr_1.4.0       pillar_1.3.1     rlang_0.3.4      rematch_1.0.1   
## [13] lazyeval_0.2.2   rstudioapi_0.10  Matrix_1.2-14    rmarkdown_1.12  
## [17] labeling_0.3     htmlwidgets_1.3  munsell_0.5.0    broom_0.5.2     
## [21] compiler_3.5.1   modelr_0.1.2     xfun_0.7         pkgconfig_2.0.2 
## [25] mgcv_1.8-24      htmltools_0.3.6  tidyselect_0.2.5 gridExtra_2.3   
## [29] crayon_1.3.4     withr_2.1.2      MASS_7.3-51.4    grid_3.5.1      
## [33] nlme_3.1-139     jsonlite_1.6     gtable_0.3.0     magrittr_1.5    
## [37] scales_1.0.0     cli_1.1.0        stringi_1.4.3    xml2_1.2.0      
## [41] generics_0.0.2   tools_3.5.1      glue_1.3.1       hms_0.4.2       
## [45] parallel_3.5.1   yaml_2.2.0       colorspace_1.3-2 cluster_2.0.7-1 
## [49] rvest_0.3.2      knitr_1.23       haven_1.1.1